function empty = checkEmptyISAnalytical(restr,irs,phi,opt)
% Function analytically determines whether the identified set is empty.
% Approach is based on the ideas in Gafarov, Meier and Montiel-Olea (2018).
% Specifically, we consider different combinations of active constraints
% and check whether there exists a vector satisfying these constraints and
% any inactive constraints. 

% This approach is only applicable when a single column of the rotation
% matrix Q is constrained by the restrictions.

% In notation of GMM18, model is 
% Y_t = c + A_1*Y_t-1 + ... + A_p*Y_t-p + B*eps_t, eps_t ~ N(0,I).
% Restrictions are on columns of B rather than on rotation matrix, Q. 
% B = A0^(-1) = Sigmatr*Q.

% Inputs:
% - restr: structure representing restrictions
% - irs: structure containing IRs and possibly long-run (cumulative) IRs
% - phi: structure containing reduced-form VAR parameters
% - opt: structure containing model information and options

Sigmatr = phi.Sigmatr;
B = phi.B; % Note B here is the (vectorised) reduced-form VAR coefficients
eqRestr = restr.eqRestr;
signRestr = restr.signRestr;
vma = irs.vma;
lrcir = irs.lrcir;
p = opt.p;
const = opt.const;
jshock = opt.jshock;
signNorm = opt.signNorm;
                
empty = 1; % Initialise indicator for empty IS - 1 = empty, 0 = not empty

N = size(Sigmatr,1); % No. of variables in VAR
mZ = size(eqRestr,1); % No. of equality restrictions
mS = size(signRestr,1) + 1; % No. of sign restrictions (incl. normalisation)

% If there are equality restrictions only and if mZ <= N-1, then the 
% identified set is nonempty by Proposition 3 (Convexity).
if (mZ <= N-1) && (mS == 1) 
    empty = 0;
    return   
end

Sigmatrprime = Sigmatr';
Sigma = Sigmatr*Sigmatrprime;
Sigmainv = Sigma\eye(N); 

if mZ > 0 % Construct matrix of equality restrictions

    Z = zeros(N,mZ);
    
    for ii = 1:mZ

        if eqRestr(ii,4) == 1 % Restriction on A0 (i.e. B^(-1))

            Z(:,ii) = Sigmainv(eqRestr(ii,2),:)';

        elseif eqRestr(ii,4) == 2 % Restriction on A0^(-1) (i.e. B)
            
            Z(eqRestr(ii,1),ii) = 1; 

        elseif eqRestr(ii,4) == 3 % Restriction on B^(-1)*A_l (A_l in GK18)

            B = reshape(B,[N*p+const,N]);
            LagCoef = B(const+(eqRestr(ii,3)-1)*N+1:const+(eqRestr(ii,3))*N,...
                eqRestr(ii,2))';
            Z(:,ii) = (LagCoef*Sigmainv)';

        elseif eqRestr(ii,4) == 4 % Restriction on LRCIR

            Z(:,ii) = lrcir(eqRestr(ii,1),:)';

        end

    end

else % Matrix of equality restrictions is empty

    Z = [];

end 

% Construct matrix representing sign restrictions (including normalisation).

if mS > 1

    S = zeros(N,mS);

    for ii = 1:mS-1

        if signRestr(ii,5) == 1 % Sign restriction on IR
        
            S(:,ii) = vma(signRestr(ii,1),:,signRestr(ii,3)+1)'*...
                signRestr(ii,4);
            
        elseif signRestr(ii,5) == 2 % Sign restriction on A0
            
           S(:,ii) = Sigmainv(signRestr(ii,2),:)'*signRestr(ii,4);
            
        end

    end

    if signNorm == 0

        S(:,ii+1) = Sigmainv(jshock,:)'; % Normalisation on A0

    else

        S(jshock,ii+1) = 1; % Normalisation on A0^(-1) (i.e. B)

    end

else

    if signNorm == 0

        S = Sigmainv(jshock,:)'; % Normalisation on A0

    else

        S = zeros(N,1);
        S(jshock) = 1; % Normalisation on A0^(-1) (i.e. B)

    end

end

% Number of active sign restrictions to consider
sTilde = min(mS,N - mZ - 1); 

% Generate vector containing column indices representing all possible 
% combinations of sTilde of the columns of S.
sComb = nchoosek(1:mS,sTilde);
rNum = size(sComb,1); % Number of possible combinations

for rr = 1:rNum % For each set of active constraints

    % Matrix containing active constraints.
    r = [Z S(:,sComb(rr,:))];
    % Matrix containing non-active sign restrictions.
    Sna = S(:,setdiff(1:end,sComb(rr,:)));

    q = null(r'); % Compute orthonormal basis for nullspace of r'.

    for jj = 1:size(q,2) % Nullity of r' may be greater than one
        
        empty = ~(all(Sna'*q(:,jj) >= 0) || all(-Sna'*q(:,jj) >= 0)); 

        if empty == 0

            return

        end

    end

end  
        
end